Phase diagrams of Zwanzig models: The effect of polydispersity 
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The first goal of this article is to study the validity of the Zwanzig model for liquid crystals to 
predict transitions to inhomogeneous phases (like smectic and columnar) and the way polydispersity 
affects these transitions. The second goal is to analyze the extension of the Zwanzig model to a 
binary mixture of rods and plates. The mixture is symmetric in that all particles have equal volume 
and length-to-breadth ratio, k. The phase diagram containing the homogeneous phases as well as 
the spinodals of the transitions to inhomogeneous phases is determined for the cases k — 5 and f 5 
in order to compare with previous results obtained in the Onsager approximation. We then study 
the effect of polydispersity on these phase diagrams, emphasizing the enhancement of the stability 
of the biaxial nematic phase it induces. 

PACS numbers: 64.70.Md,64.75.+g,61.20.Gy 



I. INTRODUCTION 

Although the Zwanzig model was introduced long ago 
as a very simplified model to study the isotropic-nematic 
transition in liquid crystals 0], the determination of its 
phase diagram, including all inhomogeneous bulk phases, 
has never been carried out so far. With the help of the 
scaled-particle theory (equivalent to a y3 expansion), de- 
veloped for a variety of hard particle fluids |2| includ- 
ing the hard parallelepipeds with restricted orientations 
which define the Zwanzig model the phase behavior 
of the homogeneous bulk phases has been investigated. 
Also the the isotropic-nematic interface has been accessi- 
ble to calculations through a smoothed-density approxi- 
mation [j| , which consists in evaluating the bulk free en- 
ergy at some smoothed density, following the same recipe 
of Tarazona for hard spheres |^ . 

The usual second virial approximation (exact for freely 
rotating models but only approximate for restricted 
orientation models) has recently been applied to the 
Zwanzig model in order to elucidate the bulk phase be- 
havior of a polydisperse mixture of hard rods Q ■ With 
an inhomogeneous version of this approximation different 
interfaces of a monodisperse hard rod fluid , as well as 
its polydisperse counterpart Q, have been studied. 

Unlike for hard objects with free orientations — like 
spherocylinders, which have been extensively simulated 
[9|, there is only one simulation of the Zwanzig model 
in a square lattice [To| . for parallelepipeds with dimen- 
sions 5xlxZ? (with D =1, rods, or 5, plates) 0. The 
main reason for this lack of simulation data is that any 
Monte Carlo movement which includes reorientation of 
parallelepipeds to any of the three directions has a very 
low acceptance ratio due to the huge overlap between 



particles as the length-to-breadth ratio and the density 



increase. 
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The second virial approximation of the free energy of 
the Zwanzig model has also been employed to investi- 
gate the phase diagram of symmetric mixtures of rods 
and plates Stimulated by theoretical calculations 

made in the early 70's which show that a binary 
mixture of rods and plates can stabilize a biaxial nematic 
phase (in which the symmetry axes of particles of differ- 
ent types point along mutually perpendicular directions) , 
van Roij and Mulder studied the relative stability of this 
phase against a nematic- nematic phase separation [Tl|. 
The authors calculated the phase diagram for different 
length-to-breadth ratios k = L/D. They showed that 
for k > 8.8 the biaxial phase is stable in a relative small 
window above a multicritical point. How the topology of 
these phase diagrams is modified when a scaled-particle 
approximation (which includes higher order virial coeffi- 
cients) is used instead of the second virial is one of the 
open questions which we will try to answer here. Another 
one is the influence of polydispersity in the phase behav- 
ior of the mono and bidisperse Zwanzig model. This is 
most relevant since with a recently derived a fundamental 
measure functional (FMF) for hard parallelepipeds 0] 
(equivalent to the scaled-particle theory for homogeneous 
phases) , we have shown that polydispersity may enhance 
the thermodynamic stability of a biaxial nematic phase 

There are recent experiments on true mixtures of hard 
rods and plates which show a very rich phase behavior, 
including triple coexistence between an isotropic and two 
nematics — each rich in one of the species — as well as in- 
homogeneous phases, like a columnar phase 0]. This 
phase behavior has been quantitatively accounted for us- 
ing the Parson approximation of the free energy func- 
tional of binary mixture of rods and plates in the Onsager 
limit 0]. Nevertheless, the long searched biaxial phase 
has not been found in these experiments. An unavoidable 
ingredient of experiments is polydispersity. In Ref. [T^j 
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we have shown that this otherwise undesirable element 
may cause the stabilization of the biaxial nematic phase. 
It is very important, though, that the system maintains 
rod-plate symmetry, in contrast with what happens in 
the existing experiments [l5j . An open question in this 
work is the topology of the phase diagram at high packing 
fractions. Addressing this problem requires to study the 
possible phase transitions to inhomogeneous phases — like 
smectic, columnar and solid phases. The calculation of 
phase equilibria in polydisperse liquid crystals when one 
of the bulk phase is inhomogeneous is thus far a theoret- 
ical challenge. There is only one such work, which makes 
use of an approximate functional for length-polydisperse 
parallel hard cylinders to make two important predic- 
tions: the existence of a terminal polydispersity beyond 
which the smectic phase is no more stable, and the en- 
hancement of the columnar phase stability These 
results have been confirmed in simulations of freely ro- 
tating hard spherocylinders in the Onsager limit |18| . 

The paper is organized as follows. In Sec. Ill Al we 
describe both the model and the fundamental measure 
functional we will use to describe its free energy (A); 
the equations for phase equilibria between phases of a 
polydisperse system (B); and the formalism to deter- 
mine spinodal instabilities in mono, bi or polydisperse 
systems (C). Section ITTT1 describes the phase diagrams of 
the Zwanzig model without (A) and with polydispersity 

(B) , as well as the bidisperse rod-plate model without 

(C) and with polydispersity (D). Finally we conclude in 
Sec. El 



II. THEORY 
A. Model 

Let us consider a length-to-breadth polydisperse mix- 
ture of uniaxial oblate and prolate parallelepipeds, with 
their symmetry axes pointing along one of the three co- 
ordinate axes. Let us fix the volume of any particle to 1; 
thus if A = L/D is the length-to- breadth ratio (with L 
the length and D the breadth) of a parallelepiped, then 



L = A 2 / 3 , D = A^ 1 / 3 



(1) 



Let us define p v {X) to be the density distribution function 
of the species parallel to the v (= x, y, z) axis, and let the 
total number density to be 



/>oo 

/ dXp(X), p(A) = 5>"(A). 
Jo 



(2) 



i labeling different components, and the excess part, 
<& oxc , can be approximated by the FMF for hard par- 
allelepipeds [jjj, namely 

$ cxc = -n ln(l-n 3 ) + ^^ + 7 ^-, (4) 

l-n 3 (l-n 3 y 

where the weighted densities {n a } have the form 

n = Pi = P, < = ]T [ P Di + p\ {Li - A)] , (5) 



n 3 



P-L^DI n v 2 = J2 \p L i ~ Pi ( L i - A)] A- (6) 



Specializing the above expressions for our polydisperse 
mixture and using Eq. Q we obtain 



/>OC 

$id = V / dXp v (X) [hip" (A) - 1] 
„ Jo 



(7) 



- p l n( l_ /3) + ^^± + _ikit_, ( 8 ) 
A ±i/3[( AT i_ 1 )^ (A)+p(A) ]. (9) 



The pressure can be obtained from its definition 

011 = V / dX b"(A)3"(A)]-#, (10) 
„ Jo 



<r(A) 



6 P »(xy 

which in this particular case yields 



(11) 



3ii - 9 + + 2 n^£+ (u \ 
p i-p + (i-P) 2 + (i- P r { } 



B. Phase equilibria between homogeneous phases 

To obtain the phase equilibria we follow the general 
procedure already reported elsewhere 0, |2l|. Suppose 
that among the n coexisting phases there are n/ isotropic, 
njv nematic and ub biaxial phases (n = nj + + ns)- 
The global density distribution (the parent distribution) 
is fixed to be 



P(A) = P h(X), / dXh(X) = 1, 



(13) 



so total mass conservation of each species is expressed by 
the lever rule 



The (temperature reduced) free energy density of a 
multicomponent mixture is given by fiFV" 1 = $ = $id+ 
<&oxc where $id is the ideal part, whose exact form is 



pU^Pi-V: 



(3) 



(14) 



with /O a (A) the total density distribution of phase a and 
7 the fraction of the total volume it occupies (there is 
the obvious constraint ^2 a ^/ a = !)• Minimizing the free 



3 



energy density <E> with respect to the fraction of particles 
with lcngth-to-breadth ratio A oriented along the v axis, 
i.e. p u a {\) = PaW/PaW, and using Eq. fHfl. the equal- 
ity of the chemical potentials of each species in different 
phases 



/3 Mo (A)=5>*(A)^(A), 



(15) 



leads to the following expressions for the coexisting den- 
sities 



p»(X)=P Q h(X) 



,-KW 



Efc 7h E T e " 



(16) 



The parent distribution we are going to use throughout 
the paper is characterized by 

h(X) = A- 1 [C/(A/«) + (1 - C)/(A/s)] , (17) 
f(z) = Kvia)- 1 exp [-(a/2) (z 2 + z' 2 )}, (18) 

where K v (a) (a > 0) is the ^th-order modified Bessel 
function, and k > 1. This choice is motivated by its 
rod-plate symmetry for £ = 1/2. The best way to ap- 
preciate this symmetry is changing the variable to In A: 
/i(lnA) = Xh(X) has two identical humps (wider the 
smaller a) centered at \nn (rods of "typical" aspect ra- 
tio k) and — In K (plates of typical aspect ratio k _1 ). 
The parameter < ( < 1 allows one to tune the over- 
all composition of the mixture, since the molar fraction 
of the rods is given by x r (Q — f. h(X)d\ [and that 
of plates by x p (() = 1 — x r (Q, of course]. Thus we 
can select polydisperse rods (plates) by setting £ — 1 
(£ = 0). The moments of this distribution are given by 
(A m ) = K m/2 (a)K (a)- 1 [C« m + (1 - C)^ m ], explicitly 
showing the symmetry of the mixture. A quantitative 
characterization of the polydispersity can be given if we 
determine the dispersion in L and D as obtained from 
h(X) for ( = or 1. This yields 




(19) 



where v = 2 for A/, and v = 1 for Aj> 

The number of independent moments in the set {p, 
which completely determine each coexisting phase is 3, 
5 and 7 for the isotropic, nematic and biaxial phases re- 
spectively (which amounts to a total of 3n/ + 5njy + 7ns 
unknowns). They can be obtained through the defini- 
tions @ and 10 , together with the distribution functions 
(|16fl . The remaining unknowns — the n — 1 independent 
7 Q 's — are calculated from the equality of pressures in ev- 
ery phase. This leaves the global dilution, Pq, as the 
control parameter. Alternatively, we can fix an external 
pressure, IIo, and eliminate Pq in terms of this new con- 
trol parameter. As Pq increases, the fractions of volume 
of each phase change; thus, for practical purposes, it is 
computationally simpler to use one of the 7 Q 's as control 
parameter and obtain Pq as a function of it. 



A particularly important case is the two phase coex- 
istence between a phase (say /3) which fills the whole 
volume (cloud phase) and an incipient new phase (say a) 
which fills an infinitcsimally amount of volume (shadow 
phase). In a cloud-shadow coexistence the parent P(X) 
coincides with the distribution function of the cloud 
phase. Then p and Ei/£± are nxe d f° r this phase, so 
the number of unknowns reduces by 3 (the number of 
constraints). 



C. Spinodal instabilities with respect to 
inhomogeneous phases 

In order to be as general as possible in developing the 
formalism let us consider the problem of finding the con- 
ditions for the stability of an arbitrary multicomponent 
system against spatial modulations of the densities pi 
(i label the species). This way we can replace at the 
end Ei — * f° r the monodisperse Zwanzig model and 
Ei - * Ei/ / ^ for the polydisperse one, to obtain the re- 
spective spinodal instabilities. The reason for so proceed- 
ing is that, as we will show later, the spinodal equations 
can be expressed in two different ways, each of which 
suits one of the two cases. 

The condition of equality of chemical potentials be- 
tween two phases is equivalent to the minimization of 
the grand potential 



n 



J dr |fc B T$({n Q }) - ^MiPi(r) j , 



(20) 



at fixed chemical potentials of each species pi, with re- 
spect to the density profile of each species Pi(r). We 
specify this grand potential for a general model whose 
free energy density $ has an excess part depending on 
the densities only through certain weighted densities 



n Q (r) 



(r), 



(21) 



as it is the case of Rosenfeld's fundamental measure the- 
ory [l^. In the definition of n a the symbol "*" stands 
for the convolution of two functions, i.e. / * g(r) = 
fdr'f(r')g(r-r'). 

We are concerned with the case when one of the coex- 
isting phases is homogeneous and the other one is inho- 
mogeneous, with its density profile consisting in a small 
perturbation around the homogeneous phase with den- 
sity distribution pi, i.e. 

Pi(r)=p i [l + e i (r)}, e< (r) < 1. (22) 

The result of minimizing l|20|) with respect to Pi(r) can 
be cast as 



pi(r) = pi exp 



E 



[A0 O 



(r) 



(23) 
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where A(j) a (r) = </> a (r) — <fi a and <fi a = d<fi/dn a . We 
are implicitly using magnitudes without spatial variable 
arguments to denote homogeneous phase quantities. Ex- 
panding A<p a (r) to first order in e,(r) yields 



A(f> a (r) = 51 Pi [ e J 



(r) 



where <j> a p — d 2 $/(dn a dnp). Inserting lt^|) in [t2"3l and 
expanding the exponential again to first order in (r) we 
arrive at the integral equation 



i( r ) = ~ K$ 53 ft [ £ i * W a 



7 /3 



(r), (25) 



<*,0 J 

which is handier to write in Fourier space 



eM) + E E P^MH(<1) = 0, (26) 

a-fl j 

where e and ft are respectively the Fourier transforms of 
e and u) [with the usual definition /(q) = J (ire 4qr /(r)]. 

Equation l|26[) gets clearer if written in matrix form. 
Defining the functions 



rn. 



,-(q) = 53^0^(q)n^(q), (27) 

Q,/3 



the matrices M(q) = (my(q)) and P = (pjSij), and the 
vector e(q) = (e^q)), Eq. (|2~6*)l becomes 



I + M(q)-P e(q) = 



(28) 



where I is the identity matrix. This linear system has 
nontrivial solutions if and only if 



D(p, q) = det I + M(q) • P 



= 0. 



(29) 



This equation is suitable for the cases in which there is 
a small number of components, like the monodisperse 
Zwanzig model with only three components: v — x, y, z. 

Alternatively we can multiply (12(:i[) by piQL(q), sum 
over i and define the functions 



%(q) = 53ft e ^ q )°7( q ) ; 

i 

then Eq. (|26[1 becomes 

Uj{q) + 51 4>apn ai {c\)up(<\) = 0, 

where we have introduced new functions 
n a <y(q) =5Z^ fi «( q ) fi 7( < l)- 



(30) 



(31) 



(32) 



Equation l|31|) can also be written in matrix form, but 
now the indices run over the set of weights, not the 
components. Hence, defining the matrices <3E> = (dap) 



and N(q) = (n Q7 (q)), and the vector u(q) = (u a (q)), 
Eq. jSj becomes 



I + N • $ u = 



(33) 



(24) Again this system has nontrivial solutions if and only if 



D{p, q) = det I + N • <f> = 



(34) 



This alternative characterization of the spinodal is suit- 
able for polydisperse systems (actually, it is a general- 
ization of the formalism developed in Ref. |2(j for ho- 
mogeneous phases). The reason is that it replaces inte- 
gral operators by finite matrices, in which the number of 
components is limited by the number of weights of the 
theory. For instance, for the polydisperse Zwanzig model 
with a free energy functional given by the FMF derived 
m Ref. 0, functions n a p take the form 



M(q)=53 / dA^(A)^(q,A)^(q,A), 



(35) 



where a, j3 = 0, (lx, ly, lz), (2x, 2y, 2z), 3 and 

Oo(qbA) = Y[w (q k A uk ), (36) 

k 

^(q,A) = JjA ukW^{q k K vk ) , (37) 
k 

^ )fe (q,A) = k vkWl {q k k vk )Sll{c±,\), (38) 

^ fe (q,A) = [A vfc «;i(g Jb A I/Jfe )]- 1 nS(q J A) (39) 

(k = x, y, z). Here q k are the components of the vector 
q, and 

A uk = \-V* + (\W-\-V*)6 vk , (40) 
w (x) = cos(x/2), (41) 
w 3 (x) = 2sin(x/2)/x, (42) 
wi(x) = w 3 (x)/w (a;). (43) 

As there are 8 different weights, D(p,q) is in this case 
the determinant of an 8 x 8 matrix. 

Either way we express the equation, (|29|l or l|34|l . the 
spinodal instability condition of a homogeneous phase 
with density distributions p 1/ (A) (or simply p v in the 
monodisperse case) with respect to the transition to an 
inhomogeneous phase occurs at the values of the to- 
tal density p and wave vector q for which the function 
D(p, q) first vanishes. More precisely, this function has 
an oscillatory behavior as a function of q, and is positive 
as long as the homogeneous phase is stable, so the spin- 
odal corresponds to the smallest density p for which the 
absolute minimum of D with respect to q equals zero. 
This amounts to finding, for a given p, the solutions to 



VD(p,q) 



0. 



(44) 



When we later consider the rod-plate mixture we will 
have to locate the nematic-biaxial phase transition. As 
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it is a continuous transition it can also be found as a 
solution to Eq. I|34|l with q = (both phases are homo- 
geneous), using 



p u {\) = P h(X), 



(45) 



£ r e-^(A) 

as the distribution functions of the nematic phases. 

III. RESULTS 
A. Monodisperse Zwanzig model 

As a first step towards the calculation of the whole 
phase diagram of the monodisperse (pure rods or p ure 
plates) Zwanzig model using the FMF of Ref. hj we 
will determine the isotropic-nematic (I-N ± ) coexistence 
curves as well as the spinodals for the nematic-smectic 
(N ± -S ± ), nematic-columnar (N ± -C ± ) and isotropic- 
plastic (orientationally disordered) solid (I-PS) transi- 
tions, as a function of the length-to-breadth ratio k. The 
+ and — superscripts label prolate (rods) and oblate 
(plates) parallelepipeds respectively. If the transitions 
are first order, the location of the coexistence curves will 
differ from that of the corresponding spinodals. However, 
as the main concern of this paper is the effect of polydis- 
persity, we defer the calculations of coexistence with and 
between inhomogeneous phases to a forthcoming publi- 
cation. 

We have obtained the spinodals for this model by solv- 
ing Eqs. (EH) and HH) using the monodisperse version of 
(|45|l . If we choose the nematic director parallel to the z 
axis, then the N ± -S ± spinodals can be calculated setting 
q = (0, 0, q), whereas the N ± -C ± ones follow from taking 
q = ((7,0,0) or q = (0, q, 0) (both are equivalent due to 
the nematic symmetry). For the I-PS spinodal all three 
previous vectors give the same result. 

Figure shows the results of these calculations. We 
observe in this figure that there is I-N^ transition only 
for k < 0.210 and k > 5.02. At those two limiting values 
of k the N + phase is destabilized by a S + at a packing 
fraction p = 0.280, and the N~ by a C~ at a packing 
fraction p = 0.261 (remember the particle volume has 
been set to 1). Also, the N + -S + spinodal is always be- 
low the N + -C + , and the N _ -C~ spinodal is always below 
the N~-S~ one. This is what intuition tells us but, as 
we will see later, it is a peculiarity of the monodisperse 
system. Despite this, all these spinodal curves converge 
asymptotically to the same packing fraction, p — 0.314, 
as k — > oo or k — > 0. This is precisely the packing frac- 
tion at which the continuous freezing transition occurs in 
a system of parallel hard cubes (k = 1) |22|. The reason 
for this is that upon increasing k the number of rods with 
orientation perpendicular to the director becomes vanish- 
ing small, and then the system is, after rescaling the z 
direction, almost equivalent to a system of parallel cubes. 
The same holds for plates upon decreasing n. What this 




FIG. 1: Phase diagram (packing fraction, p, vs. length-to- 
breadth ratio, k) of the monodisperse Zwanzig model. Solid 
and dashed lines represent the isotropic (I)-nematic (N*) co- 
existence curves. The dashed lines correspond to the values 
for which the homogeneous phases are unstable with respect 
to the inhomogeneous ones. Dotted lines represent the spin- 
odals of the different homogeneous phases: each one is labeled 
with the corresponding phase: smectic (S ), columnar (C ) 
or plastic solid (PS). 



means is that most likely the N + -S + and N _ -C~ transi- 
tions will be metastable with respect to freezing in a large 
portion of the phase diagram. For 0.182 < k < 4.93 we 
have found an I-PS spinodal instability. The peculiarity 
of this curve is that it exhibits strong oscillations as the 
aspect ratio changes. These oscillations reflect the pack- 
ing efficiency of randomly oriented parallelepipeds as a 
function of their size (the better the packing the lower 
the curve). 

The available simulation results for freely rotating hard 
spherocylinders show that the I-S + and N + -S + transi- 
tions begin at K = 4.1 and 4.5 respectively (notice 
that for hard spherocylinders the length-to-breadth ratio 
is k = £ + p ). On the other hand, simulations of hard cut 
spheres show that for k = 0.2 there is an I-C~ transition 
(the isotropic phase might instead be a peculiar 'cubatic' 
phase) and for k — 0.1 a N~-C~ one We can see that 
despite the different particle geometry and the restricted 
orientations of the Zwanzig model, the agreement with 
the threshold k's at which spatial instabilities destabi- 
lize the homogeneous phases is rather good. Also the 
qualitative picture is similar: elongated rods form smec- 
tics, while flat disks form columnars, and more symmet- 
ric particles form solids instead of smectics or columnar. 
Thus, this simple model seems to capture the essence 
of the entropy driven phase transitions between phases 
with different symmetries and their relation to particle 
anisotropy. 

The only existing simulation of the Zwanzig model 
has been performed on a lattice model of parallelepipeds 
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with length-to-breadth ratios 15 : 3 (re — 5) and 3 : 15 0.5 
(re = 0.2) [ljj. F° r these re's the authors find I-C~ for the 
discs and and I-DS + for the rods, where DS + stands for 
a novel phase called 'discotic smectic' by the authors (in 
this phase the axes of the particles point perpendicular q ^ 

to the normal of the smectic layers while there is no ori- 
entational order within the layer) . The packing fractions 
they obtain for the transitions are higher than those of 
the corresponding spinodals of Fig^ but they should de- 
crease upon decreasing the lattice spacing, as the results 0.3 
for the freezing of parallel hard cubes on a lattice (oc- 
curring at p = 0.568 for edge-length 2 lattice spacings, 
at p = 0.402, for edge-length 6 lattice spacings, and at 
p = 0.314 for the continuum) illustrate j2J]. All this is q ^ 

again in agreement with the phase diagram of Fig. ^ 1 2 1 1 10° 1 1 10 
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B. Polydisperse Zwanzig model: Unimodal 
distribution 

Let us study now the effect polydispersity has on the 
phase behavior of the Zwanzig model discussed above. 
We have introduced polydispersity through the unimodal 
parent probability density resulting when setting £ = 1 
(rod-like parallelepipeds) or £ — (plate-like paral- 
lelepipeds) in function (|T7|) . Before reporting the results 
we have obtained, a few words on polydisperse phase di- 
agram plots are on purpose. 

Generally speaking, plotting polydisperse phase dia- 
grams would require an infinite dimensional space, for as 
we vary the fraction of total volume occupied by the coex- 
isting phases they change their composition. A phase di- 
agram like that of Fig^plots densities vs. aspect ratios, 
thus carrying no information whatsoever about composi- 
tions. So we must restrict ourselves to a given composi- 
tion, and this is the parent one, since all pure phases have 
this composition. This means that only cloud lines can 
be plotted, delimiting regions where we can only know 
that the system is decomposed into two or more coex- 
isting phases. The result resembles the coexisting lines 
of monodisperse phase diagram, but the meaning is com- 
pletely different. Although we could connect cloud points 
in the two lines delimiting the coexistence region, they 
are definitely no coexisting states. 

Similar considerations hold for spinodal lines. We can 
represent the spinodal of a system in a pure phase having 
the parent distribution as composition, but if the system 
reaches a cloud point an incipient shadow phase with 
a totally different composition — hence a totally differ- 
ent spinodal line — coexists with it. If the density of this 
shadow phase is above its own spinodal line, the phase 
transition will not be stable. 

Having all this in mind we have solved the coexistence 
equations of the isotropic and nematic phases and deter- 
mined the spinodal lines of inhomogencous instabilities 
for a system with the parent distribution l|17|) . The re- 
sults for the choice a = 1 (corresponding to polydisper- 
sity parameters = 0.288 and Ad = 0.143) are shown 



FIG. 2: Phase diagram for polydisperse Zwanzig model. The 
distribution of aspect ratios is unimodal with the maximum 
at re. Polydispersity is Al = 0.288 in length and Ad = 0.143 
in breadth. The lines have the same meanings as in Fig. 
Full circles mark the points above which the nematic cloud 
phase is unstable against columnar ordering; empty circles 
mark the points above which the isotropic cloud phase coex- 
ists with a nematic shadow phase unstable against columnar 
modulations (the spinodal line of this nematic shadow, com- 
puted with the corresponding aspect ratio distribution, is not 
plotted in the figure). 



in Fig.lU 

Apart from the usual broadening of the isotropic- 
nematic transition already observed for this model in 
Ref. H we find a clear enhancement of the stability of 
the homogeneous phases i.e. the spinodal instabilities oc- 
cur at higher densities. This effect is more pronounced 
for the N ± -S ± or N ± -C ± transitions. Another remark- 
able feature of Fig. [21 is that the columnar phases are 
by far more stable than the smectic ones for any aspect 
ratio. 

This general behavior is in agreement with the re- 
sults obtained from both density functional theory for 
length-polydisperse parallel hard cylinders |17| and com- 
puter simulations for length-polydisperse freely rotat- 
ing hard spherocylinders in the Onsager limit [lg. Al- 
though in Ref. [17] some strong approximations were 
made, such as the decoupling of the density profile as 
p{r,L) = g(L)p(r), and the effect of fractionation be- 
tween the nematic and columnar phases was ignored, 
the results agree qualitatively with the simulations: both 
show a terminal polydispersity beyond which the smec- 
tic phase is no more stable, being replaced by a columnar 
phase (which tolerates a higher degree of polydispersity) . 
The nematic-columnar transition is first order. While 
fractionation is indeed negligible in the nematic-smectic 
transition, the smectic-columnar transition clearly segre- 
gates long rods into the columnar phase. 

It is interesting to notice that in the limits re — > oo 
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or none of the spinodal lines tend to the density of the 
freezing of parallel hard cubes (n = 1). The reason is that 
no trivial scaling can be applied to a polydisperse system 
of perfectly aligned parallelepipeds to transform it into 
a system of parallel hard cubes. On the other hand the 
S + and S - spinodals tend to the same limiting density 
when k — > 0,oo, and so do the C + and C~ spinodals. 
The second virial term of the free energy is dominant for 
large (small) k's; this and the rod-plate symmetry of the 
chosen parent distributions explains this behavior. 



C. Binary mixture of Zwanzig rods and plates 

We have calculated the phase diagram (coexistence of 
homogeneous phases as well as their spinodal instabili- 
ties to inhomogeneous phases) for the binary mixture of 
rods and plates with n — 5 and 15. The results for k = 5 
are shown in Fig. [3J Apart from the usual I-N ± tran- 
sitions and N~-N + phase separations which coalesce in 
a multicritical point we see that the nematic phases in- 
crease their stability with respect to the inhomogeneous 
ones as the relative composition of the fluid tends to the 
equimolar mixture. Another important feature is that 
the biaxial nematic is not stable. 

All these characteristics are also present in the phase 
diagram obtained using a second virial approximation 
[Tl|. but there are important differences as well. The 
most prominent one is that the whole phase diagram is 
shifted down in density, with all transitions and spinodals 
appearing below p — 0.4 (those obtained with a second 
virial approximation occur above that value This 
result should not surprise if one takes into account that 
in a y3 expansion the excess free energy is written in 
terms of the variable y = p/(l — p), which grows faster 
than p. Another striking difference with respect to the 
second virial approach is the rod-plate asymmetry of the 
phase diagram, a consequence of the inclusion of higher 
virial coefficients — which do not share the symmetry of 
the second one. 

We pass now to describe the loss of stability of the 
nematics with respect to inhomogeneous phases. The in- 
tersection between the N + line of the I-N + coexistence 
and the N + -C + spinodal in the rod-rich part of the phase 
diagram (marked in the figure with a full circle; the open 
circle corresponds to its coexisting isotropic phase) is con- 
sistent with the existence of a first order phase transition 
between the I phase and an inhomogeneous (presumably 
columnar) phase. (Recall that simulations on the lattice 
find a first order I-DS + transition for parallelepipeds with 
n = 5 10]). Nevertheless a definitive answer about the 
relative stability between these inhomogeneous phases — 
including also the solid — can only be given after carrying 
out coexistence calculations. 

Increasing the fraction of plates (which in the C + phase 
have their principal axes oriented perpendicular to the 
director) seems to favor the smectic alignment of rods. 
This results in the intersection of the N + -C + and the 
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FIG. 3: Phase diagram of a binary mixture of rods and plates 
with the fraction of rods and p coincides with the 

packing fraction. Spinodal and coexistence curves are repre- 
sented by dotted and solid lines, respectively. Dashed lines 
are coexistence curves in the unstable region. Fhe crossings 
between the coexistence curves and the spinodals are repre- 
sented with full circles, and with empty circles their coexisting 
phases. Full squares mark the crossover between columnar 
and smectic phases along the spinodals. 



N + -S + spinodals (marked with a full square on the right 
of Fig.[3J). In the second virial approach the N + -S + spin- 
odal is always below the N + -C + one. At the plate-rich 
part of the phase diagram the nematic exhibits a spin- 
odal to the C~ phase, as usual in discotic fluids In- 
creasing the fraction of rods beyond a threshold value 
(indicated with a full square on the left part of the phase 
diagram) makes the smectic more stable than the colum- 
nar, which agrees qualitatively with the results from the 
second virial approach [ll"|. 

The phase diagram for n — 15 appears in Fig. 0] Its 
main difference with the previous one is the presence of a 
thcrmodynamically stable biaxial nematic phase in an in- 
verted triangular right above the multicritical point. The 
window is limited to the left and to the right by continu- 
ous phase transitions to uniaxial nematic phases, and to 
the top by a B-N + coexistence which becomes a N~-N + 
phase separation as density increases. As discussed in 
Ref. |ll| . the driving mechanism which determines the 
preference of this system for phase separation instead of 
biaxial ordering is the larger exclude volume between un- 
like particles compared to that of like particles (the rod- 
plate excluded volume divided by the rod-rod one scales 
as k 2 / 3 for large k). When the gain in free volume com- 
pensates the loss in mixing entropy (and this strongly 
depends on concentration and composition) phase sepa- 
ration occurs. 

The slight asymmetry of the phase diagram is again 
due to the presence of higher virial terms in the free en- 
ergy. The second virial approximation predicts a small 
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FIG. 4: Phase diagram of binary mixture of rods and plates 
with k — 15. Lines and symbols have the same meanings as in 
Fig. 03 Notice that given the continuous nature of the B-N* 
transitions, the spinodals are the true transition lines. 



region of triple coexistence between the N _ , N + and B 
phases right on top of the biaxial phase region. The ab- 
sence of this triple coexistence region in Fig. [|] is due to 
the asymmetry of the phase diagram. Simulations of a 
mixture of prolate and oblate ellipsoids |25| show a simi- 
lar asymmetric phase diagram, only tilted in the opposite 
direction (there is a B-N~ demixing on top of the biax- 
ial nematic, instead of the B-N + of the present model), 
probably due to the different particle geometry. Despite 
the asymmetry, we should admit that for k = 15 there 
is good overall agreement between the second virial ap- 
proach and ours. 

We have not studied the changes an asymmetry in the 
mixture (a difference in the volume of rods and plates, 
for instance) would induce in the phase diagram. Intu- 
itively, asymmetry acts against the stability of the biax- 
ial phase, something that is confirmed by a recent work 
|26j which studies a binary rod-plate mixture of freely 
rotating oblate and prolate cylinders in the Onsager ap- 
proximation. This work explicitly shows that the region 
of stability of the biaxial phase decreases with increasing 
the asymmetry of the mixture (although, according to 
this work, this phase can be found in highly asymmetry 
mixtures). 

With respect to the inhomogeneous phases, Fig. [I] 
shows a remarkable behavior: although pure hard rods in 
the N + phase undergo a spinodal instability to a smectic 
one, by adding a tiny fraction of plates the instability 
takes place to a columnar phase. This phenomenon has 
already been observed in a different system, namely, in 
simulations of binary mixtures of parallel hard sphcro- 
cylinders j2]|- Although pure parallel hard spherocylin- 
ders exhibit a continuous nematic-smectic transition, a 
binary mixture of them with length-to-breadth ratios 2 



and 2.9 show a first order nematic-columnar transition 
instead. This result was explained by the poorer pack- 
ing of rods of different length in the smectic phase as 
compared to that of rods of the same length. Borrowing 
the argument, in our rod-plate model the plates can fit 
into the interlayer spacing with their axes parallel to the 
smectic director as long as there are few of them; but 
upon increasing their molar fraction some of them are 
forced to get into the smectic layers with their axes per- 
pendicular to the smectic director, thereby destabilizing 
the smectic phase. 

Finally, at the other extreme of phase diagram the N~ 
phase losses its stability to a columnar phase of plates at 
densities lower than those of the N + -C + spinodal. 

D. Polydisperse Zwanzig model: Bimodal 
distribution 

In this section we are going to study the effect of poly- 
dispersity on the phase behavior of the previous rod-plate 
mixture. In Ref. [T^j we have shown how polydispersity 
can stabilize the biaxial phase even for mixtures with 
relative small aspect ratios — like k = 5, for which Fig. EH 
shows that the biaxial phase is absent. 

Apart from the increase in mixing entropy that poly- 
dispersity carries, the other mechanism behind the en- 
hancement of stability of the biaxial phase is the decrease 
of the ratio of the average exclude volumes of like and un- 
like particles. The excluded volume between two parallel 
rods with lengths Li and breadth Di (i = 1, 2) is 

v rr = (L 1 +L 2 ){D 1 +D 2 f , (46) 

and the excluded volume between a rod with length L\ 
and breadth D\ and a plate with length L 2 and breadth 
D 2 , with axes perpendicular to each other, is 

v rp = {L 1 +D 2 ){L 2 + D l ){D 1 + D 2 ). (47) 

According to (JIJ, if the volume of all particles is set to 
one, in terms of the length-to-breadth ratios A^, Li = 

2/3 1/3 

A i / and Di — X i . Taking a double average of v rr 
over J], [h{\i)Q{\ - 1)] and of v rp over h{\i)] 0(Ai- 
1)6(1 — A2) [h(X) is the parent distribution (jT7Jl ]. and 
fixing £ = 1/2 (the equimolar composition) we arrive at 

((v rp )) \ +m 2 _ 1 + m\ + 2(m_im 2 + mim_ 2 ) , . 

77 \T ~ i ' \ ") 

{{Vrr)) 5 + 4m_imi + 2m^ 2 m 2 

where mp = / 1 °°dA/ l (A)A' 3 / 3 . It can be shown numeri- 
cally that, for a given k, the ratio l|48|l decreases with 
polydispersity (i.e. with decreasing a). Analytic ex- 
pressions for this ratio can be obtained in the limit of 
high particle anisotropy («> 1) and high polydispersity 
(a <C 1) with the constraint an 2 <C 1 (in terms of the 
parameters A„, v = L,D, this constraint is equivalent 
to 1 <C A„ <C \/ln «, implying that the fraction of cubic- 
like particles is vanishing small). Using the asymptotic 
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FIG. 5: Phase diagram of a polydisperse mixture of rods 
and plates with k = 5 and length and breath polydispersi- 
ties A_t=0.610 and Ad ~ 0.302. Part (b) is a zoom of the 
region right above the biaxial phase. 



oo 0.000 0.000 6.9 

2.00 0.216 0.107 6.6 

1.00 0.288 0.143 6.3 

0.50 0.374 0.185 5.9 

0.25 0.470 0.233 5.1 

0.10 0.610 0.302 4.0 



TABLE I: Threshold values of the length-to-breadth ratio, k* , 
for the appearance of a biaxial nematic phase as a function 
of polydispersity, expressed in terms of the parameters a, Al 
and An- 



expressions 



771/3 



K /3/3 



*/V 6 (-lna)' 



we obtain 



(Kp)) 

({v rr }) 



(-lna) 1/2 . 



v = L,D 



(49) 



(50) 



with CV a positive constant. This asymptotic relationship 
explicitly shows both the scali ng o f this ratio discussed 
for the pure rod-plate mixture [lj and the exponential 
attenuation of this scaling with polydispersity. 

We have estimated the threshold value of k (denoted 
k*) beyond which the biaxial phase begins to be thermo- 
dynamically stable, as a function of polydispersity. The 



results are summarized in Table [I] Without polydisper- 
sity this value is k* — 6.9, smaller than the second virial 
estimation ^l| n* — 8.8. This is an indication that three 
(and higher) body correlations increase the stability of 
the biaxial phase. As we can see in Table U n* decreases 
upon increasing polydispersity, an illustration of the en- 
hancement of stability of the biaxial ordering induced by 
polydispersity. 

In Fig. [S] we plot the phase diagram of this system for 
k = 5 and a = 0.1 (corresponding to = 0.610 and 
Ad = 0.302). According to Table [I] for this parame- 
ter choice the biaxial phase is thermodynamically stable, 
and this is clearly shown in Fig. |5Ja), where we can see 
a small window of biaxial nematic just above the multi- 
critical point. Figure Efb) is a zoom of the upper border 
of this window. There is a narrow region of N~-B co- 
existence limited from above by a N - -N + coexistence. 
The dividing line (the dotted line in the diagram) rep- 
resents a continuous transition between the B and the 
N + phases which coexist with the N - . This line has no 
analog in the binary rod-plate mixture. To obtain it one 
has to solve the equations for two phase coexistence at 
all values of 7 in 116|) (the volume occupied by one of the 
coexisting phases) and then find the one for which the 
B-N + transition occurs. 

The high density part of the phase diagram of Fig.JSJa) 
was not determined in Ref. 14]. The most remarkable 
feature it shows is that pure phases are hardly stable: the 
diagram is dominated by coexistence regions. The reason 
is that polydispersity is so high that the mixture has, 
beside rods and plates, a significant amount of cubic-like 
particles. This favors entropic phase separation in the 
three different phases I, and N + . 

In the upper part of Fig. |Sf a) there is a region of triple 
coexistence I-N~-N + . The lines limiting this region were 
obtained by solving the coexistence equations for the case 
when one of the coexisting phases is a shadow phase (i.e. 
it occupies a vanishing fraction of the total volume). As 
an illustrative example, if we take 71 = 0, 7n- = 7 and 
7n- =1—7 and solve the coexistence equations for the 
N = 15 unknowns (see the discussion about phase equi- 
libria in polydisperse systems in Section III B|) we obtain 
the curve joining the two points marked with full squares 
in Fig.[Sfa). The other two curves limiting the triple co- 
existence region were calculated selecting 7 N + = (the 
left one) and 7 N - = (the right one). 

We have checked that the nematics below the triple 
region are always stable with respect to spatial density 
modulations of any kind. In order to estimate the pack- 
ing fraction at which these instabilities occur we solved 
Eqs. (|34|l and 144|) using the density distributions p v (X) 
of all the coexisting phases along the borders of the triple 
coexistence region. The left and right open circles corre- 
spond respectively to the packing fractions at which the 
coexisting N + and N~ (the shadow phase in both cases) 
lose their stability to the C + and C~ phases. Above 
these points the phase diagram will possibly include co- 
existences between four phases: I, N - , N + and C~ or C + ) 



10 



0.4 



0.3 - 



(a) 




FIG. 6: Distribution functions at the N^-I-N + (a) and the 
Nj^-I-N - (b) triple points of the phase diagram of Fig. |5] 
(marked with full squares). The subscript 'sh' denotes shadow 
phases. Solid lines represent the shadow phases, dotted lines 
the I phase and dashed lines the N + (a) and N~ (b) phases. 



(the experiments of Ref. 15] find a similar scenario). 

Figurc[S]shows the density distribution functions of the 
phases coexisting at the points marked with full squares 
in Fig. EJa). As already discussed, the curves illustrate 
the high fraction of cubic-like particles partly responsible 
for the strong demixing this system exhibits. The figure 
also illustrates the strong fractionation which takes place 
between the coexisting phases: the isotropic phase is rich 
in cubic-like particles, the N~ in plates and the N + in 
rods. 

In contrast, the phase diagram of the polydisperse rod- 
plate mixture for k — 15 does not change qualitatively 
compared to the binary mixture (cf. Figs. [7| and the 
size and shape of the biaxial regions are similar and the 
general topology of the phase diagram is basically the 
same. The only important difference (minor in terms of 
the size of the portion of phase diagram it involves) ap- 
pears just above the biaxial: one can observe both B-N + 
and B-N~ coexistence (the former occupying a larger 
region), as well as a triple coexistence zone separating 
them [see Fig. Efc)]. As for n — 5, the B-N* coexis- 
tences become N + -N~ coexistence through a second or- 
der phase transition of the B phase to the N ± [dotted 
lines in Fig. [7Jb) and (c)]. 

The N~ and N + lose their stability to columnar phases 
at relative high packing fractions [see the dotted lines in 
the upper part of Fig. 0a)]. In both cases it is the N~ 
phase which first becomes unstable with respect to C~ at 
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FIG. 7: Phase diagram and two details of a polydisperse rod- 
plate mixture with re = 15 and length and breadth polydis- 
persity Al = 0.288 and Ad = 0.143, respectively. 



the N + -N~ coexistence (the N~ cloud on the left, marked 
with a full circle in Fig. [7Ja), and the N~ shadow on the 
right, marked with an empty circle in this figure). 

Figure[S]shows the density distribution functions corre- 
sponding to the point marked by a full square in Fig. Etc). 
The distributions of the N + and B phases are very simi- 
lar, except for the fact that the N + has a slightly higher 
proportion of rods (and correspondingly less of plates) 
than the B phase. It is clear from the figure that the 
proportion of cubic- like particles is negligible, so this mix- 
ture can be regarded as a true (polydisperse) rod-plate 
mixture. 



IV. CONCLUSIONS 

We have studied the effect of polydispersity on the 
phase diagrams of several variants of the Zwanzig model 
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FIG. 8: Density distribution functions at the point marked 
with a full square in Fig. Jf \c) ■ The solid, dashed and dotted 
lines correspond to the cloud B, shadow N + and shadow N~ 
phases, respectively. 

for liquid crystals. We have first determined the spinodal 
instabilities of the homogeneous phases to the smectic, 
columnar or plastic solid phases for one component flu- 
ids of rods or plates. There is a general qualitative agree- 
ment with phase diagrams of oblate and prolate particles 
obtained from Monte Carlo simulations 0, H] ■ We have 
also shown that even a small degree of polydispersity 
has dramatic effects on the transitions to highly ordered 
phases (like the smectic or the columnar). In particular, 



we have shown that upon increasing polydispersity the 
transition to the columnar phase of rods preempts the 
N + -S + transition, in agreement with simulations |23j . 

Next we have obtained the phase diagrams of sym- 
metric rod-plate binary mixtures with length-to-breadth 
ratios k — 5 and 15. While for large k's the differences 
between these phase diagrams and those calculated with 
a second virial approximation |llj| are only quantitative, 
for small k's both methods differ qualitatively. For in- 
stance, all the transitions occur at packing fractions well 
below those obtained with the second virial approach, 
and sometimes the relative stability between different in- 
homogeneous phases changes. 

We have shown that the introduction of polydisper- 
sity in a rod-plate mixture in a very symmetric way en- 
hances the stability of the biaxial phase even for as small 
an aspect ratio as k = 5. For this case the amount of 
polydispersity needed to stabilize the biaxial phase is so 
high that the mixture becomes very unstable with re- 
spect to phase separation at practically any composition. 
There are large regions of three phase and possibly four 
phase coexistence (I-N"-N+ and I-N"-^-^). These 
results agree with what is observed in experiments with 
rod-plate mixtures with a high degree of polydispersity 

m 
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